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Abstract 

We develop a new powerful method to reproduce in silico single-molecule manipulation experiments. We demonstrate that 
flexible polymers such as DNA can be simulated using rigid body dynamics thanks to an original implementation of 
Langevin dynamics in an open source library called Open Dynamics Engine. We moreover implement a global thermostat 
which accelerates the simulation sampling by two orders of magnitude. We reproduce force-extension as well as rotation- 
extension curves of reference experimental studies. Finally, we extend the model to simulations where the control 
parameter is no longer the torsional strain but instead the torque, and predict the expected behavior for this case which is 
particularly challenging theoretically and experimentally. 

Citation: Carrivain P, Barbi M, Victor J-M (2014) In Silico Single-Molecule Manipulation of DNA with Rigid Body Dynamics. PLoS Comput Biol 10(2): el 003456. 
doi:l 0.1 371/journal.pcbi.l 003456 

Editor: Shi-Jie Chen, University of Missouri, United States of America 

Received August 13, 2013; Accepted December 11, 2013; Published February 20, 2014 

Copyright: © 2014 Carrivain et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Funding: This work was funded by the Institut National du Cancer, PL8I0 program, grant INCa_5960 and the Institut National de la Sante et de la Recherche 
Medicale, grant MICROMEGAS PC201104. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the 
manuscript. 

Competing interests: The authors have declared that no competing interests exist. 
* E-mail: victor@lptmc.jussieu.fr 



This i.s a PLOS Computational Biology Methods article. 

Introduction 

The mechanical and topological properties of DNA and 
protein-DNA assemblies are of primary importance in many 
biological processes, including transcription, replication, chroma- 
tin organization and remodeling. Since techniques have become 
available enabling the manipulation of single-molecules [1,2], a 
large amount of experimental data have been accumulated on the 
mechanical response of DNA and protein-DNA assemblies under 
stretching forces and twisting torsions, in particular from optical 
and magnetic tweezers experiments [2-5]. In magnetic tweezers 
experiments, a DNA molecule is grafted at one end to a coverslip 
and at the other end to a magnetic bead. The bead is trapped in 
the magnetic field of a pair of magnets that may be translated, thus 
exerting a varying force on the bead. Moreover the pair of 
magnets may be rotated at a certain number of turns, thus 
constraining the finking number of the DNA molecule. After the 
stretching force and the number of turns are applied to the bead, 
the only physical variable that can be directiy measured is the 
DNA extension, i.e. the distance between its two ends. Therefore 
the interpretation of the experimental results requires an 
important modeling effort, particularly in the more complex cases 
where DNA is associated with proteins, as for instance in 
chromatin assemblies [6,7]. Although theoretical approaches 
may be successful in some cases [8-10], simulations are often 
crucial tests of the proposed model vafidity, when they are not the 
unique possible way of dealing with the system complexity. 



In this spirit, we aim to develop an efficient tool to manipulate 
single-molecules in silico reproducing optical and magnetic 
tweezers experiments. This task is challenging since the DNA 
model should have precise specifications to reproduce the behavior 
of DNA accurately. We need to: (i) model a polymer, i.e. an 
articulated chain; (ii) reproduce the effective diameter of DNA 
(depending on electrostatic conditions) and, when proteins are 
present, have the possibility to model their shape and steric 
hindrance; (Hi) deal with collisions, especially in order to 
reproduce DNA supercoUed structures (plectonemes) and steric 
effects in DNA-protein assemblies; (iv) reproduce DNA twisting 
and bending elasticities; (v) include statistical mechanics features to 
account for temperature and thermal motion. Beside these 
essential points, we also wish to simulate the system dynamics, 
which may be important in some cases, e.g. when hysteresis is 
observed under magnetic tweezers [7] or for in vivo chromosome 
dynamics experiments in the cell nucleus [11,12]. 

This ambitious list of specifications is beyond the reach of 
traditional simulation approaches where particles interact through 
2-body potentials (as in Molecular Dynamics or Monte Carlo 
simulations [13,14] with a given force field). The need to deal with 
frozen degrees of freedom in coarse grained modeling may be 
addressed through holonomic constraints, as in the SHAKE 
algorithm [15,16], where an iterative approach is adopted. 
However, collision detection and steric hindrance may only be 
accounted for in this scheme by introducing additional steps. More 
recently, non iterative algorithms have been developed [17], that 
subsequently led to the development of new powerful tools, called 
"physics engines". These have been designed by the engineering 
and robotics communities to reproduce the dynamic behaviour of 
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Author Summary 

Video game techniques are designed to simulate rigid 
body dynamics of macroscopic bodies, e.g. characters or 
vehicles, in a realistic manner. However they are not able 
to deal with temperature effects, hence they are not able 
to deal with molecules. In order to extend these powerful 
techniques to molecular modeling, we implement here 
Langevin Dynamics In an open source library called Open 
Dynamics Engine. Moreover we add a "global thermostat" 
to this Langevin Dynamics, which accelerates the simula- 
tion sampling by two orders of magnitude. With these 
radically new simulation techniques, we prove that we can 
accurately reproduce single-molecule manipulation exper- 
iments in silico, in particular force-extension as well as 
rotation-extension curves of reference experimental stud- 
ies. The method developed here represents an unparal- 
leled tool for the study of more complex single molecule 
manipulation experiments, notably when DNA interacts 
with proteins. Furthermore the simulation technique that 
we propose here has all the functionalities required to 
tackle the nuclear organization of chromosomes at every 
length scale, from DNA to whole nuclei. 

articulated systems of rigid bodies. Physics engines are acquiring 
an increasing importance, notably in the fields of computer 
graphics and video games, where they are now widely used to 
simulate rigid body motion under realistic conditions and in real- 
time. Open Dynamics Engine (ODE) is one of the most popular 
rigid-body dynamics open source library for robotics simulation 
applications [18]. As other physics engines, ODE simulates the 
kinematics of articulated systems by using permanent joints that 
impose holonomic constraints, instead of bond potentials. The 
same method is used to manage collisions: when overlapping 
between bodies is detected, a temporary joint is locally created that 
reproduces the action of the contact forces, without the need for 
explicit permanent 2-body interaction potentials (see section 
"Materials and Methods" for details on how ODE manages joints 
and collisions). 

These extremely efficient simulators haven't, up to now, been 
used in statistical mechanics. Although well adapted to mechanical 
simulations, physics engines lack coupling to a thermal bath. The 
main novelty of our approach is the implementation of Langevin- 
Euler equation in the ODE software. Moreover we improve the 
simulation efficiency of this Langevin dynamics by extending the 
"global thermostat" algorithm designed by Bussi and ParineUo in 
2008 [19] to physics engines. This algorithm allows much faster 
yet unbiased sampling of the phase-space. As a first step toward 
simulating DNA-protein assemblies, we focus here on bare DNA 
and show how to perform in silico single molecule manipulation of 
DNA. 

Materials and Methods 

Introduction to physics engines 

In rigid body dynamics simulations run with ODE, the state of a 
system consisting of N rigid bodies is described by the positions r, 
of their centres of mass, a quaternion representation of their 
orientations y,-, and their linear and angular velocities v,- and co,- 
respectively. These velocities are collected in the column vector 
V = {vi,(Ui, . . . VN,(Off}^ . We use the superscript T to denote the 
transpose of a vector or a matrix everywhere in this article. The 
vector C = MV then collects all linear and angular momenta, 
where M = {inil^i . . . mffl,Iff} is a 6Nx6N block diagonal 



matrix whose elements are the mass matrices in, l and inertia 
matrices X, of the N bodies (with 1 the 3x3 identity matrix). The 
Newtonian dynamics equation then reads C = J- where the 
generalized force is a vector collecting forces and torques 
applied to the system. These forces and torques may be external, 
due for example to gravity or magnetic fields, or internal, as a 
consequence of the mechanical constraints between the rigid 
bodies that make up the system. 

Most notably, in articulated systems, as is the case of polymers, 
rigid bodies are connected by mechanical joints. A joint is a 
relationship that is enforced between two bodies so that they can 
have only certain positions and orientations relative to each other, 
and ODE provides different types of joints according to the kind of 
articulation that has to be implemented, e.g. ball-and-socket, 
hinge, slider or universal. 

Mathematically a joint imposes some holonomic constraint 
between both connected bodies. Such a constraint is an equation 
that reads 3 = 0 where d is the distance between both joint 
bearings, e.g. the center of the ball of one body and the the center 
of the socket of the other one. The constrained distance 3 is purely 
geometrical, depending only on the relative position and 
orientation of both jointed bodies. The position and orientation 
of each of the N bodies the articulated system is composed of 
depend on time /. Therefore the constraints ^ = 0 of the articulated 
system can be derived with respect to time to get the kinematic 
constraints in the form JV = 0 where we introduce the jacobian 
matrix of constraints J (see subsection "Exact solution for X when 
there are no collisions" for a detailed example). This velocity- 
based description is used in ODE as in most game/physics 
engines. 

So, mechanical joints exert reaction forces and torques on the 
joint bearings. These internal mechanical constraints can be 
collected into a generalized constraint force J-c which, by virtue of 
the principle of virtual work T^V = 0, reads T c = J^'^ where 1 is 
a vector of Lagrange multipliers that precisely accounts for the 
reaction forces and torques coming from the joint bearings [16]. 
The Newtonian dynamics equation therefore reads C, = J-g + Tc 
where and Tc stand for the external and internal contributions 
to the generalized force respectively. As the constraint force reads 
Tc = J^'^, the Newtonian dynamics equation becomes an 
equation for X in the form: C = J-e-\- J^^'k. 

Solving this equation for "k should moreover satisfy the 
holonomic constraints ^ = 0 at every timestep However the 
discretization used in the numerical calculation results in errors on 
3 so that 3(t) is generally not equal to 0. Then, in order to have 
<5(?-|-A?) = 0 at the next timestep, the kinematic constraint J^V 
should be adapted accordingly. Indeed 3{t -\- t!^t) = 3{t) + 
^(?)V(?-f A?)A/ according to the Euler semi-implicit integration 
scheme which is used in velocity-based algorithms. Hence 

J'(OV(?-f Ar)= — But then this imphes that the kinematic 

constraint is not equal to zero at time J-l-A?, i.e. 
^(?-|-Af)V(?-(- A?)#0, so that the joint bearings will continue to 
move apart afterwards. In order to keep both 3 and J^V close to 
zero at every timestep, ODE introduces an error reduction 
parameter kerp in the kinematic equation i7(?)V(? + A/) = 

— k„p-^ [18]. This parameter has to be adjusted to some 

optimal value between 0 (no correction at all) and 1 (complete 
correction of (5 in one timestep). However setting k„p = 1 is not 
recommended since, as said above, this would imply that the joint 
bearings will continue to move apart afterwards with maximal 
velocities. ODE recommends values between 0.1 and 0.8. 
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In addition to ke,p, ODE introduces a second ingredient to 
soften tlie rigid constraints by allowing the violation of the 
constraint equation by an amount proportional to the restoring 
force "k. More explicitly, a "constraint force mixing" diagonal 
matrix /Ccfm = ^f/m^ is defined, such that i7(0V(? + A?) = 



kcfm'kit) (implicit integration) [18]. This is equivalent 

k,- 

to introducing a spring-damper system (spring constant of 



and damping constant of - 



,\-ke, 



k. 



in, 



the joint bearings; this can be understood as analogous to a bead- 
spring model. Nevertheless there is a major difference between 

this effective spring and a regular spring: the term — k^rp 

constrains the velocity whereas a regular spring constrains the 
acceleration. As a result, no energy is stored in this effective spring, 
at odds with a regular spring which stores an averaged energy kgT 
(see below subsection "Preliminary tests of validity and perfor- 
mance of the global thermostat" along with the histograms of 
energy in Figure 1). 

In particular, ODE uses a powerful software called libccd [20] to 
detect collisions between two convex shapes. Whenever overlap- 
ping is detected between two rigid bodies, ODE attaches a 
temporary joint between them called a "contact joint". Defining 
vector C\ (resp. Ci) that connects the center of mass of body 1 (resp. 
2) to the contact point and denoting n the common normal to both 
bodies at the contact point (directed from 2 to 1), the kinematic 
constraint imposed by the contact joint would read 
/i-'^(vi+£Bi®ci)— n-^(v2 + 'y2®C2)=0 in the perfect case when 
the holonomic constraint imposed by the contact joint reads 
exactly <5 = 0. However, in practice S is not equal to 0 because of 
discretization errors, hence the kinematic constraint imposed by 
the contact joint actually reads: 

J{t)V{t^At) = 

k (1) 
- ^ (n^(ri +ci -r2 ~C2)) (?) -k,,,,, (n^X) (/) 



with 



J= «^ -n 



X^n>0 



(2) 



(3) 



The right hand side of Eq.(l) deals with the already existing 
overlapping of the two bodies in contact at time t when collision is 
first detected, or with their residual overlapping while the contact 
joint exists. 

By inserting the constraint force Tc = J^^k into the equation of 
motion C=T and taking the first-order discretisation of this equa- 
tion, one can easily get the following expression to be solved for X: 



Ar 



(4) 



This equation is of the form .AX = B. Importantly, the addition 
of the term — to each diagonal term of the matrix J M^^ 



global thermostat velocity based 
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Figure 1 . Boltzmann statistics test. The distribution of l<inetic energy 
P(E) of a DNA molecule of length L = 1 jj.m at thermal equilibrium is 
plotted as a function of the dimensionless kinetic energy \iE for both 
global (left panel) and local (right panel) thermostats. The DNA molecule 
is composed of A' = 300 rigid cylinders of radius r = 2 nm and / = 3.34 mn. 
The other parameters of the simulation are given in Table 1. According to 
Eq. (33) the number of degrees of freedom is dof = (sN — n^ with n,. = 2>N 
the number of non-redundant holonomic constraints (3 per ball-and 
socket joint, hence dof = 900). The factor T{n — 1) ensures the 
normalization of P(E) (F is the Euler Gamma function). 
doi:10.1371/journal.pcbi.1003456.g001 



provides a symmetric positive definite matrix A, thus gready 
increasing the solution accuracy of Eq.(4). From this equation, the 
vector X of Lagrange multipliers, hence the constraint force J-c, 
can be determined. Then the motion solver (semi-implicit Euler 
integrator) gives the new positions and orientations of the 
articulated bodies at time t + At. It is advantageous to choose 
the Exponential Map parametrization [21] for the quaternion 
integration. 
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Solving constraints 

In general Eq.(4) has to be solved numerically and ODE has two 
algorithms to do so, one based on the Successive-Over-Relaxation 

(SOR) method [22] and the other based on the Linear 
Complementary Problem (LCP) [23]. LCP time complexity is of 
order rrv' and space complexity (memory) of order np' where m is 
the number of constraint rows [18]; whereas SOR time complexity 
is of order mNgoR where NgoR is the number of successive-over- 
relaxation and space complexity of order m [18]. Both algorithms 
have equivalent performances when m = \/NsoR- But in general 
LCP is more accurate, although much more time consuming, than 
SOR. We compared these two algorithms for a chain of length 
= 300 without noticing significant differences in the errors on S 
(error on the colocalization of joint bearings). In order to save 
computational time, we preferentially run the SOR method with a 
value of C050J? = l-7 for the relaxation factor and NsoR = WO. 
These values are different from the default values in ODE and 
work well for a linear chain of rigid bodies connected with ball- 
and-socket joints. But in some cases, the SOR method does not 
converge and we then switch to the LCP method, which always 
comerges. However in the case when there are no collisions 
between the rigid bodies the articulated system is composed of, we 
were able to derive an exact solution for k (see next subsection). 
Therefore we solve Eq.(4) according to the following scheme: 

1. if there are no collisions at time t, we use the exact solution for 

2. if some collisions are detected, then we run the SOR method 
and check for the accuracy of the solution. More specifically 

k 1 

the solution is accepted if 6+ JV(t + b.i)\l 



Ktk, 



cfm ^cfin 



ll>^ll<10-1; 

3. if \\l+-r^S+- — JV{t + At)\\/\\X\\>lO-'^ then the 



Atk, 



cfm 



^cfm 



simulation step is restarted with the LCP method. 



Exact solution for 1 when there are no collisions 

For the sake of clarity, let us first consider the example of four 
rigid cylinders of length / each connected with ball-and-socket 

joints at the extremities ± - with the first one anchored to 

some fixed point, taken as the origin of the coordinates. The vector 
ti is the tangent to the cylinder /. The jacobian matrix J 
associated with this system is tridiagonal when there are no 
collisions, in which case it reads: 
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(5) 



For each cylinder, the antisymmetric matrix f® is associated 
with the cross product ?(S) and has the property <®^= — <®: 



The transpose Jacobian matrix and mass matrix M are 
given by: 



J' 



M- 



/-I 




1 


0 


0 


\ 






/ ^ 

f0 
2 ' 




2'' 


0 


0 








0 




-1 


1 


0 








0 




2^2 


2^2 


0 








0 




0 


-1 


1 








0 




0 


^ ,® 

L 


' t® 

L 








0 




0 


0 


-1 








n 

V " 




n 


0 


2'^ 


/ 






/ mi 


0 


0 


0 


0 


0 


0 


o\ 


0 


/l 


0 


0 


0 


0 


0 


0 


0 


0 


mi 


0 


0 


0 


0 


0 


0 


0 


0 


ll 


0 


0 


0 


0 


0 


0 


0 


0 


mi 


0 


0 


0 


0 


0 


0 


0 


0 


/, 


0 


0 


0 


0 


0 


0 


0 


0 


ma, 


0 


V 0 


0 


0 


0 


0 


0 


0 


u) 



Then we get: 
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m 



-1 



-1,® 
I 

1 

-1,® 
0 



2^' 



m- 



-1 



2^2 '2 



-m 



-1 



2 3 '3 
0 



2^3 h 



-1 



0 



-m 



2 ' 4 ) 



(7) 



(8) 



and we deduce the final result for the symmetric matrix 



H = 

J--Tf -L-Tf + -L-Tf 
0 



V 



0 



(9) 



0 



ty \ 

0 / 



(6) where we define the matrix T? 

in the associated principal axis body frame as T? 



, ^f/r'rf. We then write Tf 
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(10) 



matrix J M. ^ -\ — with the following properties: 



W„ = — -Tf + 

m\ A/ 



(13) 



The vector X collects the Lagrange multipliers associated to 

each joint respectively 'k = \\^^^j^. The equation 'HX = X gives us 
a system of coupled equations on X.,-. Note that ODE solves in one 
time all the constraints of this articulated system. This is not the 
case with the SHAKE algorithm where an unconstrained step is 
first performed, before correcting the positions and orientations 
iteratively to get the constraints satisfied eventually. The term 
JM^^Te from Eq.(4) is given by: 
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(11) 



H«=-^ Tf_i + l-Tf + %^ i>\ (14) 



Hff+i = Tf and Hii+x=Hi+u (15) 

mi 



Hij = Q \i-i\>l 



(16) 



Using the following decomposition CDC^ for the matrix H 
where £ is a block lower matrix with block identity matrix on the 
diagonal and where 2? is the block diagonal matrix it is easy to 
show that Cij = 0 for |/ — !| > 1. From these we get the following 
equations: 



Vi = Hu-Ci+uViCl^y i>l 



(17) 
(18) 
(19) 



In order to solve the linear system of equations TIX = X we 
define X' ='DC^X and solve the problem CX' = X in an iterative 
way: 



where we write Fx = (r^m)m+ (r^ntjm and Fn = (T^ijt. We 
can then write the constraint forces and torques J-c = J^X 
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(12) 



X," = Xi — Cii-\ Xi-\ i>\ 



(20) 
(21) 



and we get the final solution for X by solving the problem 
X'=VC^X: 



Xf^ — Vff Xff' 



Xi = Vr%,-Ci+uh+\ i<N 



(22) 
(23) 



The method explained here is the exact solution of the problem 
TiX = X where no collisions are present in the system. With this 
exact resolution the simulation is faster than the SOR algorithm 
{NsoR = 100 and cosoR = 1 -7) with a gain of 4. 



We can now generalise the previous example to the case of a 
linear chain of rigid cylinders connected with ball-and-socket 
joints with the first one anchored to the ground. We denote H the 



Simulating DNA molecules under magnetic tweezers 

To model a DNA molecule, we build a linear chain of rigid 
cylinders of length / = 3.34nm each, corresponding to 10 base 
pairs (10 bp), which amounts to the double helix pitch. We 
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connect the cylinders to each other by ball-and-socket joints. The 
radius of the cylinders is set according to the salt buffer 
concentration of the experimental data we compare with. Indeed, 
since the DNA molecule is highly negatively charged, DNA-DNA 
electrostatic repulsion affects the double helix response in single 
molecule experiments [9,24-27]. This effect can be easily and 
implicidy included in simulations and theoretical models by 
introducing an effective DNA radius r = rQ + where tq is the 
crystallographic radius of the DNA double helix and r^/ accounts 
for the DNA-DNA electrostatic reptolsion [28-31]. It turns out that 
Vei may be set equal to the Debye length of the salt buffer 
solution. As /,rf = 10/y^ with c the salt concentration given in 
mmol and in nm, we set the effective radius to r = 2 nm in 
100 mmol monovalent salt buffer for comparison with the 
reference experimental data of Mosconi et al [32]. Alternatively 
we use an effective radius r = 4 nm to fit the experimental data 
obtained in 10 mmol monovalent salt buffer by Smith et al [1]. 

We performed all our in silico single molecule experiments with a 
DNA molecule of contour length L=l |im. The corresponding 
number of DNA cylinders in the chain is therefore N = 300. The 
DNA molecule is anchored, at one end, to a planar surface 
(mimicking the microscope coverslip), and at the other end, to a 
rotatable bead (mimicking the magnetic bead). We set the bead 
radius to L/2ti in order to prevent the DNA from looping around 
it. At both ends of the DNA chain, the rigid cylinders are tangent 
to their attachment surface. 

The final problem that remains to be addressed is how to obtain 
a correct definition of the bending and twisting behavior of DNA. 
We have solved this problem by a special choice of the connecting 
joints and by introducing appropriate restoring torques reacting to 
the bending and twisting deformations. This has been done based 
on the bending and twisting energies that are defined according to 
the usual expressions PEi,=gi,{\— COS 6) and fiE,=g,^^ /2 
respectively. The rigidity constants gh and gt are related to the 
bending and twisting persistence lengths p and t respectively, 
through the following equations: 



c(gb) = icose-) 



<COS0> = 



2p-l 
Ip + l 



(24) 



(25) 



gt = t/l 



(26) 



where C is the Langevin function (see supplementary Text SI). 
The bending angle 9e[0,7t[ and twisting angle (j>e[ — n,Ti\ are 
related to the standard Euler transformation ZXZ and are given 
by 



COS 0: 



(27) 



(28) 



(1 +f^/2)sin(^ = »ij»ii — mJiMi (29) 

where t, is the tangent vector of cylinder i, »i, a vector normal to ti 
and mi = ti®mi. These three vectors are the principal axis of 
cylinder i. 

We finally get the following expression for the global restoring 
torque between two connected DNA segments (see supplementary 
Text SI for the complete derivation of this equation): 



b+t- 



--ght\®t2- 



g4 



1 +cos( 



(30) 



We recall that, for DNA, estimates of the bending persistence 
length give /i = 50nm for 10 — 500 mmol salt buffer (see Refs. 
[1,33-35]); whereas estimates of the twisting persistence length 
give ? = 95 nm for 10— 100 mmol salt buffer [9,32]. According to 
the size / = 3.34nm of the unit cylinder we find gb = l5-5 and 
^/ = 28.4. 

Langevin dynamics and global thermostat 
implementation 

Although well adapted to mechanical simulations, ODE lacks 
coupling to a thermal bath. As physics engines impose to deal with 
dynamics equations including inertial terms, in particular for 
computing constraint forces (collected in Tc), we need to turn to 
some implementation of stochastic isothermal molecular dynamics 



10 n 




I I I I I I I 

1 10' 10^ 10^ 10* 10^ io' 



Figure 2. Global thermostat efficiency test. Complementary autocorrelation function 1 — Crr(E*t) of the end-to-end distance of a DNA 
molecule simulated with the global Langevin thermostat (black) compared with the same function simulated with the local Langevin thermostat 
(red). Z*T is the dimensionless lag-time with S* the thermostat coupling frequency. The DNA molecule is composed of A' = 300 rigid cylinders of 
radius r = 2 nm and / = 3.34 nm. The other parameters of the simulation are given in Table 1. Cylinders are connected by ball-and-socket joints. The 
chain ends are free to diffuse here. 
doi:1 0.1 371/journal.pcbi.1 003456.g002 
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Table 1. Main set of parameters used for the DNA model and for the numerical simulations. 





Entity 


Parameter 


Typical value 


Definition 


DNA 

parameters 


r 


2nm 


DNA cylinder effective 
radius in 100 mmol salt 






4 nm 


DNA cylinder effective 
radius in 10 mmol salt 




/ 


3.34 nm 


DNA cylinders length 
(10 basepairs} 




N 


300 


number of DNA cylinders 


(corresp. to a contour 
length L = 1 (im) 




P 


50 nm 


bending persistence length 
in 10-100 mmol salt 




t 


95 nm 


twisting persistence length 
in 10-100 mmol salt 




gb 


15.5 


bending rigidity constant 


(corresp. to bending persistence length : p = 50 nm) 




gt 


28.4 


twisting rigidity constant 


(corresp. to twisting persistence length : t = 95 nm) 




mo 


1.162 10-^^ kg 


mass of each DNA cylinder 
(10 basepairs) 


Time units 


to 


1.7660 I0-'° s 


natural system time unit to = l^molkBT 


Simulation parameters 


At 


0.000592 to 


simulation time step 




Z* 


10 


thermostat coupling frequency 






0.8 


ODE error reduction parameter 






lO-'ro/mo 


ODE constraint force mixing parameter (hard) 






1.7 


SOR relaxation factor 




^SOR 


100 


number of successive-over-relaxations 



doi:l 0.1 371 /journal.pcbi.l 003456.t001 



in order to thermalize the system: isothermal to simulate the system 
at constant temperature, stochastic to ensure ergodicity. The 
corresponding algorithms are all related to Langevin dynamics 
and can be cast into local and global thermostats. In local 
thermostats, such as standard Langevin dynamics, a correction 
force including both a frictional term and a stochastic term is 
exerted on each particle to drive the system to the canonical 
distribution at a prescribed temperature. Global schemes of 
Langevin dynamics are designed to minimize the perturbation 
introduced by the thermostat on the Hamiltonian trajectory (so 
called "disturbance" as defined originally in [36]), hence on the 
dynamical properties, such as autocorrelation functions, and 
related quantities, such as difiFusion coefTicients. In these globally 
applied thermostats the stochastic term of the correction force 
acting on each particle is proportional to the momentum of that 
particle. Two main global algorithms have been designed so far: (a) 
Stochastic Velocity Rescaling methods, most notably the "global 
thermostat" introduced by Bussi and Parrinello [19], (b) the Nose- 
Hoover Langevin thermostat [37]. Here we first show how to 
implement Langevin-Euler equation in the ODE software. More- 
over we show that the global thermostat introduced by Bussi and 
Parrinello is so remarkably adapted to this implementation that it 
improves quite significantly the sampling efficiency with respect to 
local Langevin dynamics (by two orders of magnitude in typical 
situations), while preserving the time-dependent properties such as 
autocorrelation functions. The sampling efficiency is defmed as 
usual as the number of independent configurations generated 
during the time necessary to reach thermal equilibrium. 



To begin with, we add to the "mechanical" forces T an 
additional, thermal contribution Q = — T.C + SW containing a 
frictional term — E£ and a random force vector EW. E is the 
matrix of the coupling frequencies to the thermostat, 3 the matrix 
of white noise amplitudes and W a generalized vector of 
normalized and independent Wiener processes. 3 and E are 
related through the fluctuation-dissipation theorem, which reads 
here 



(31) 



where /? = 1 2" with T the temperature of the thermal bath and 
where the superscript ★ denotes that the matrices E and 2 are 
chosen to be diagonal in the principal axis body frame (where the 
matrix M is diagonal by definition). For simplicity, we choose to 
fix all the S* to a common frequency S*. Note that (S*)^' is the 
relaxation time of the thermostat, i.e. the autocorrelation time of 
the kinetic energy (see supplementary Figure SI). 

We then improve the sampling efficiency of this Langevin 
dynamics by extending the "global thermostat" algorithm 
designed by Bussi and ParineUo in 2008 [10] to physics engines. 
This algorithm allows faster yet correct sampling of the phase 
space in the canonical ensemble. However, it is designed for the 
translational degrees of freedom only. In order to apply it to an 
articulated rigid body system, we therefore have to extend it to the 
rotational degrees of freedom and adapt it to the ODE software. 
To this aim, we replace the traditional Langevin-Euler correction 
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force (local thermostat) 5 by a corresponding global version Q, 
which reads 



(32) 



Eq.(32) shows that Q is proportional to jC, so that the stochastic 
force and torque globally associated with the thermostat is in the 
same direction as £.. Hence, a free particle, i.e. a particle not 
connected to any other particle, will move on a straight line 
between two collisions. Note that, nevertheless, the particle will 
undergo Brownian motion along diis straight trajectory. The 
global version of the Langevin dynamics minimizes the distur- 
bance induced by the thermostat on the HamUtonian trajectory 
(equal to M^^Q according to its definition in Ref. [19], but 
extended here to the rotational degrees of freedom), nevertheless 
retaining the same thermalization speed as usual Langevin 
dynamics (see supplementary Figure SI). 

When used in the framework of a velocity-based algorithm such 
as ODE, the global thermostat presents a remarkable advantage. 
This is because, in this case the global Langevin contribution Q is 
decoupled from the constraint forces, in the sense that it cancels 
out in the equations for J^^- More precisely, with our definition of 
Q (see Eq.(32)), the contribution JM^^Q to the term JM^^T in 
Eq.(4) is always zero. In other words, G not only minimizes the 



100.0 -I 
50.0 



10.0 
5.0 



1.0 
0.5 - 



A numerical WLC 

analytical WLC fit to exp (1 0 mM salt) 

o LpNA = 300x10 bp 10 mM salt 

analytical FJC 

+ Ldna= 10x300 bp 




<z„>/L 

Figure 3. Comparison to force-extension curves. Red circles: 
dimensionless stretching force Pfp as a function of the mean relative 
extension £ = <zo>/i- Here the radius of the cylinders is i' = 4nm 
corresponding to the 10 mmol monovalent salt buffer used in [1]. The 
other parameters of the simulation are given in Table 1. For 
comparison, the black solid line reproduces the analytical Worm-Like- 
Chain force-extension approximation formula [35]. Black triangles 
correspond to a numerical fit of the exact Worm-Like-Chain model 
[34] with the same persistence length p = 50 nm. Blue crosses: we also 
show simulations in the limit case l = 2p, with no torsional rigidity 
(g/=0) and no collisions, and compare it to the theoretical force- 
extension curve of a Freely-Jointed-Chain (FJC, blue solid line). The 
statistical error bars on the simulation points are all smaller than the 
symbol size. 

doi:1 0.1 371/journal.pcbi.1 003456.g003 



disturbance of the HamUtonian trajectory M^^G, but also does 
not disturb the generalized constraint force J-c- Both effects 
cooperate to achieve a dramatic acceleration of the simulation 
sampling, that is, in the case of our model, approximately 100 
times faster than with the local thermostat (see below "Preliminary 
tests of validity and performance of the global thermostat" and 
Figure 2). Importantly this acceleration is compatible with the 
correct computation of dynamical properties, such as autocorre- 
lation functions. 

Parameter settings of our implementation of ODE 

In all DNA simulations presented in this article, we choose the 
length of the cylinders as the unit length /q = /, the mass of the 
cylinders as the unit mass OTq = m and the unit of thermal agitation 
ksT as the unit of energy Eo=kBT, from which we deduce the 
unit of time to = /q ^/mo/Eo. The complete set of parameters of 
our simulations is given in Table 1. We also choose to deal 
coUisions with a restitution coefficient equal to 1 without surface 
friction. Hence, when two rigid bodies collide, the constraint force 
X imposed by the contact joint that temporarily connects them is 
directed along their common normal n at the contact point. We 
finally choose an error reduction parameter A:„.^ = 0.8 and a 

_9 ^ 
mo' 



constraint force mixing parameter kcfk-^ 



10- 



Results/Discussion 

Preliminary tests of validity and performance of the 
global thermostat 

We start to validate our methodology by simulating a DNA 
molecule without any constraint applied on the bead (neither 
stretching nor twisting). To this aim, we first check the 
equipartition theorem. When the system is at thermal equilibrium, 
its temperature is related to the kinetic energy through the 
equation 



2<£> = <£^M-'£> = {6N-nc)kBT 



(33) 



where ric is the number of non-redundant holonomic constraints. 
This relation is standard since 6N — tic is just the number of 
degrees of freedom (dof) of the system. Moreover the distribution 
of the kinetic energy of the system at thermal equilibrium follows a 
Boltzmann law and therefore reads: 



P(E)- 



T(n-\) 



(34) 



with n = dof /2— 1. We checked this relation for a DNA molecule 
of length L=\ |j.m coupled to the two different Langevin 
thermostats, local and global respectively. The resulting histograms 
are shown in Figure 1, confirming that: (i) the kinetic energy is 
correctly sampled at thermal equilibrium with both thermostats, (ii) 
there is indeed no energy stored in the joints, although these have 
been softened by effective springs (see above the error reduction 
parameter kerp in subsection "Introduction to physics engines"). 

We then quantified the simulation sampling efficiency by means 
of the autocorrelation function Crr(t) = <(/f(r-|-T)/f(r))/(2L/i) of 
the end-to-end distance R of a DNA molecule of length L= \ jim 
coupled to the two different Langevin thermostats, local and global 
respectively. Here the average is performed over the time t and T 
denotes the lag-time of the autocorrelation function. A demon- 
stration of the performance of the global thermostat in terms of 
relaxation rapidity is given in Figure 2. Fitting the exponential 
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Figure 4. Comparison to experimental extension-rotation curves. IVIean relative extension £ = <z^>/_L as a function of the imposed overtwist 
a for different stretching forces / in the range of 0.35— 1.80 pN. We superimpose on to the simulation results (symbols) the experimental results from 
[32] (lines). The statistical error bars on the simulation points are smaller than the symbol size. 
doi:1 0.1 371 /journal.pcbi.1 003456.g004 



decrease of both relaxation curves shows that with use of the global 
thermostat we reach the saturation value at S*T =10* whereas this 
value is reached at S*t= 10* in the case of the local thermostat, 
thus resulting in an acceleration factor of about 100 for this system 
composed of N = 300 articulated rigid bodies. Note that the 
dimensionless lag-time S*T is equal to (S*?o)(A?/?o)(t/A?) ~ 
0.006 ftstep (see Table 1), with risiep — z/At the corresponding 
number of time steps. Then a typical run using the global thermostat 
is of the order of tens of millions of time steps, whereas it is of the 
order of billions of time steps with usual Langevin dynamics. A 
striking illustration of the sampling acceleration provided by the 
global thermostat is also given in supplementary Video S 1 . 

We also compute the tangent-tangent correlation function ('f fy) 
along the polymer with both local and global thermostats. No 
significant deviations were formd between both thermostats. Results 
obtained with the global thermostat are plotted in supplementary 
Figure S2 along with the corresponding theoretical curves. 

A simple calculation shows that the tangent-tangent correlation 
function decreases as exp( — \j — i\l /p), from which one can 
calculate the average bending <cos 9). With the DNA persistence 
length /i = 50nm, this quantity amounts to 0.9355, to be 
compared to the result from a fit of the simulation curves, giving 
<^cos 0)j„„ = 0.9359. The same comparison can be done for the 
twist angle (with / = 95nm), for which the simulation average 
'v'^^Xim = 0.0352 matches the theoretical value 0.0352. These 
comparisons show that our simulation results are in very good 



agreement with the analytical formulae, thus validating (i) our 
implementation of the bending and twisting rigidities and (ii) the 
correct sampling of the DNA conformation space by means of the 
global Langevin thermostat. 

Comparison with the experimental DNA stretching 
response 

We then simulate reference force-extension curves, both 
theoretical and experimental. We thus perform simulations at 
given stretching force fe along the z-axis (normal to the DNA 
anchor surface), and without torsional constraints on the magnetic 
bead. In order to fit the experimental data obtained by Smith et al 
[1] in 10 monovalent salt buffer, we set here the DNA radius (i.e. 
the radius of the unit cylinders) to r = 4 nm. The resulting force- 
extension curve is given in Figure 3 where we plot (red circles) the 
dimensionless stretching force jifp as a function of the dimension- 
less mean relative extension s,= (z(j} j L. Here <Zo) denotes the 
mean DNA extension, i.e. the mean distance between the bead 
and the anchor surface, at zero torsional constraint. For 
comparison, we also plot (black solid line) the analytical Worm- 
Like-Chain (WLC) interpolation fitting curve proposed by 
Bouchiat et al [35] as well as the numerical solution of the 
WLC model (black triangles) obtained by Marko and Siggia with 
the same persistence length p = 50 nm [34] . The simulation 
reproduces pretty well the WLC behavior, thus validating our 
implementation of the DNA bending rigidity. Note that, at low 
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25 -1 o f=0.74 pN/r=2.00 nm 
A f=0.91 pN/r=2.00 nm 
+ f=1.13 pN/r=2.00 nm 



20 - 



15 - 



10 - 



5 - 



.1 



4 

4 









± 



4 







_Q o.^ L. 



0 



0.00 



0.05 



0.10 



0.15 



Figure 5. Torque computation. Twisting torque F as a function of the average overtwist <(t> for stretching forces / = 0.74, 0.91 and 1.13 pN. 
Symbols are the simulation results. Horizontal lines show the values of the critical torques estimated from the experimental results at the 
corresponding stretching forces [32]. Oblique black dotted line: experimental results in the pure extended state (low <(t». Oblique blue dotted line: 
asymptotic behavior of the twisting torque as a function of the simulated average overtwist in the pure plectonemic state (high <(t»: 
^«^T»=^'B^(27r//;o)T(<(T>-o■p) with T = 49.60 + 0.77 nm and (7^-0.0610 + 0.001; /io = 3.57 nm is the pitch of the DNA double helix. 
doi:1 0.1 371 /journal.pcbi.1 003456.g005 



forces, the extension saturates at a value greater than zero becau.se 
of the impenetrable ground and magnetic bead that both confine 
the DNA molecule. This effect is more pronounced than in the 
experimental curve [1] because the ratio of the bead radius to the 
DNA length is higher in our simulations. 

In Figure 3, we also show the results obtained in the limit case 
l = 2p, for which gA=0 (see Eqs.(24— 25)), and when there are no 
collisions. In this case we expect to observe a Freely-Jointed-Chain 
response (with iV=IO segments of 300 bp each). The analytical 
force-extension relation for FJC is given by the well-known 
expression as a Langevin function C{2Pfp) and it is also 
reproduced in Figure 3. Again, the simulation results are in very 
good agreement with the theoretical formula. 

Comparison with the experimental DNA torsional 
response 

More interestingly, magnetic tweezers also allow the application 
of a torsional strain on a single DNA molecule at constant 
stretching force. This torsional strain is equal to the number of 
turns of the magnetic bead around the z-axis due to the rotation of 
the magnets. The number of turns of the bead is also equal to 
ALk, the variation of the linking number of the DNA double helix 
with respect to the intrinsic twist of the DNA double helix 
Lko = Z//Ao with Ao = 3.57nm the pitch of the DNA. And we 
define as usual the DNA relative overtwist as f7 = ALk/Lko. 



Simulations at constant strain (fixed overtwist). We 

simulate the rotation-extension behaviour of our DNA model by 
imposing the bead rotation and compare to experimental data 
from [32], where the mean relative extension £=<(z^)/L is given 
for different stretching forces as a function of the fixed relative 
overtwist a (Figure 4). An example of the simulated dynamics at a 
stretching force / = 0.74pN and ALk =15 is shown in supple- 
mentary Video S2. Again excellent agreement is observed between 
the experimental bell-shape curves (also called "hat curves") and 
the simulated curves. We went further to check the validity of the 
DNA radius which is set here to r = 2 nm according to the 1 00 
monovalent salt buffer used in these experiments. To this aim we 
refer to a series of papers by Neukirch and co-workers [9,24,25] 
where they showed that the linear part of the "hat curves" can be 
expressed as a function of the supercoiling radius p, the 
superhelical angle 6 and the ratio j of the twisting and bending 
persistence lengths, as 

<2o> r ^cos20 J Hsinie"' ^ ' 

where p and 9 depend on the force applied on the bead. 

The comparison between the experimental estimations of the 
three parameters p, 6, and 4np/{H sinlO) deduced from [32] and 
the corresponding results from our simulations are given as 
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Figure 6. Buckling behavior. Bucl<ling instability with stretching force of/ =0,74 pN and six different torques 9.0,9.5,9.6,9.7,9.8,9.9 pN nm. (Blue) 
DNA relative extension. (Red) DNA overtwist. DNA relative extension and overtwist are monitored as a function of the number of simulation steps. All 
these recordings have been obtained after thermal equilibrium has been reached. 
doi:1 0.1 371 /journal.pcbi.1 003456.g006 



functions of the applied force in the three supplementary Figures 
S3, S4, S5. The good agreement observed validates the value 
applied to the effective DNA radius r = 2 nm in 100 mmol 
monovalent salt buffer. 

Simulations at constant stress (fixed torque). In the 
characteristic rotation-extension "hat curves" the DNA extension 
decreases linearly as a function of ff once it reaches a critical point 
where the resulting applied torque crosses a critical buckling value 
Ff. This linear decrease corresponds to the formation of 
plectonemes: any additional turn beyond the buckling transition 
is absorbed into the growing plectoneme without changes in the 
torque [10]. As standard magnetic tweezers cannot measure 
torques, the buckling torque can be only indirectly deduced by 
integrating the change of the molecule extension with respect to 
the applied force [32,38]. However, new set ups have recently 
been proposed that enforce a torque and allow its measurement: 
one of them uses an angular optical trap [39-41], another one a 
magnetic nanorod coupled to a magnetic bead [42] and a third 
one a soft magnetic tweezer [43]. These innovative experiments 
confirm that the torque stays constant during the plectoneme 
formation, and allow the investigation of the dependence of on 
the applied stretching force. 

Our technique also allows us to simulate the DNA response 
when both the stretching force / and the twisting torque F are 
imposed, while the number of turns of the magnetic bead Lk is 
free to evolve. In this case we compute the average overtwist (cr) 



for frxed values of the twisting torque F. Figure 5 shows F(<fT» 
for three different values of the external force. A clear transition 
from "pure extended" DNA (left part of the curves) to "pure 
plectonemic" DNA (right part) is observed at (almost) constant 
critical torque F^, In the pure extended state, simulations and 
experimental results [32] are in very good agreement. In this 
regime F({f7)) is linear and the corresponding slope agrees with 
theoretical predictions [10]. The estimations of the critical torque 
obtained in [32] by integrating the change in the molecule 
extension are reported as horizontal lines in Figure 5. Our 
simulation correctly reproduces these values, including the 
dependence of F^. on the stretching force / throughout the whole 
range of salt concentrations explored (data not shown). 

Note that we sporacUcaUy obtained two plectonemes in the same 
DNA molecule during some simulation runs. Moreover we also 
observed that plectonemes diffuse along the DNA molecule when 
the global thermostat is started up (see supplementary Video SI), 
Interestingly, both features - multiplectonemes and plectoneme 
diffusion - have been recentiy observed experimentally [44] and 
theoretically explained [45], Note that nevertheless these obser- 
vations have been obtained with a 21-kb DNA molecule, hence 
with a much longer molecule than in our simulations (3-kb), This 
may explain why we did not observe the "multiplectoneme phase" 
described in ([45]). 

The increase of the average overtwist (ff) observed at the 
critical torque in Figure 5 corresponds to the formation of 
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plectonemes, until the entire DNA molecule is supercoiled. We 
explored this critical regime by recording, at given / and F, the 
time evolution of both the molecule length and overtwist. When 

the critical torque is approached, a clear buckling instability 
appears, with DNA fluctuating between the pure extended state, 
characterized by a smaU supercoiling and a large extension, and a 
pure plectonemic state, with opposite characteristics [10] (Figure 
6). Note the anticorrelation between supercoiling and extension. 

Beyond the buckling torque, DNA is in the pure plectonemic 
state, and an increasing torque wiU tighten the molecule 
supercoiling. Interestingly, while experimental results [32] do not 
provide measurements of the average overtwist •((t) in this pure 
plectonemic state, our simulations aUow to explore this regime. An 
example of the simulated dynamics at a stretching force 
/ = 0.74pN and twisting torque r=15pN-nm is shown in 
suppk^mcntary Video S3. We observe that our simulation results 
deviate from Marko's theoretical prediction according to which 
r(<0'» =A:5rfpiecto(27t//!o)<(T>, with fpiecto the twist persistence 
length of the plectoneme and hi) the pitch of the DNA double helix 
[10]. Instead we get an affine law which may be cast in a similar 
form: r«(7» =A:i;r(27r//io)T(<(7>-CT^) with T = 49.60 + 
0.77 nm and (7^,^0.0610 + 0.001. The presence of the reference 
value Op seems rather natural since this is the overtwist above 
which the system enters the pure plectonemic state. The meaning 
of the T length is an open question, but we notice that its value is 
close to the DNA bending persistence length p. 

Concluding remarks 

The method developed here can conveniently describe all 
physiological situations involving DNA positive supercoiling, 
which are of main importance for DNA transcription or 
replication. One limitation of the modeling developed so far is 
that extensive modifications of the double helix structure are not 
accounted ioi, e.g. base pair opening that occurs when negative 
supercoiling is applied (DNA denaturation), or the S-DNA 
transition observed at extremely high force, or the P-DNA 
transition under very positive torque. Nevertheless this method- 
ology gives an unparalleled opportunity to study more complex 
biological .systems, such as protein-DNA complexes: in particular, 
we are currently addressing the modeUng of magnetic tweezers 
response of chromatin fibres [7], by associating our in silico DNA 
model with a rigid body representation of the histone octamer. 
More recentiy, in vivo experiments also enable the measurement of 
the dynamics of chromosome bci in the cell nucleus [11,12,46]. 

Supporting Information 

Figure SI Autocorrelation function of the kinetic energy 
CeeO^*!) (symbols) of a polymer of size L = l Jim and persistence 
length p = 50 nm and for three different cylinder lengths: / = 3.34 
(circles), / = 6.68 (triangles) and /= 10.02 nm (crosses). We 
compare the theoretical exponential model exp{ — \J — i\l/p) for 
the tangent-tangent correlation to our simulation results (blue 
lines) . 
(TIF) 

Figure S2 Average tangent-tangent correlation (tjtj') (symbols) 
for a polymer of size L=\ (tm and persistence length p = 5Q nm 
and for three dilferent cylinder lengths: / = 3.34 (circles), / = 6.68 
(triangles) and / = 10.02 nm (crosses). We compare the theoretical 
exponential model exp( — [/ — /|///>) for the tangent-tangent 
correlation to our simulation results (blue lines). 
(TIFF) 



Figure S3 Supercoiling radius p estimation from simulation results 
(black circles) and from experimental results [32] (red triangles). 

(TIFF) 

Figure 84 Helical angle estimation 0 from simulation results 
(black circles) and from experimental results [32] (red triangles). 

(TIFF) 

Figure S5 Slope estimation 4np/{H sinlff) from simulation 
results (black circles) and from experimental results [32] (red triangles). 
(TIFF) 

Text SI In this appendix we address the problem of how to 

obtain a correct and computationally efficient definition of the 
bending and twisting rigidities for a polymer in general and a 
DNA molecule in particular. 

(PDF) 

Video SI Illustration of the sampling efliciency of global versus 
local thermostat. This 1 min 23 s video records the simulated 
dynamics of a DNA molecule manipulated by magnetic tweezers 
after thermal equilibrium has been reached. The stretching force is 
/ = 0.74pN and the number of turns of the bead is fixed to 
ALk= 15. From 0 to 20 s the simulation is run with Langevin local 
thermostat. Then the global thermostat is started up at 21 s and 
run untU 41s, strikingly accelerating the dynamics of the DNA 
molecule. Moreover the plectoneme starts to diffuse along the 
DNA molecule. Then the local thermostat is run again from 42 s 
to 1 min 02 s, resulting in a dramatic slowing down of the 
dynamics. And finally the global thermostat is started up again at 
1 min 03 S and run until the end of the video at 1 min 23 S. The 
DNA molecule is composed of A = 300 rigid cylinders of radius 
r = 2 nm and / = 3.34 nm. The other parameters of the simulation 
are given in Table 1. This video is also available at http:/ /vimeo. 
com/51918121. 
(AVI) 

Video S2 Simulated dynamics of a DNA molecule manipulated 
by magnetic tweezers with a stretching force / = 0.74 pN and fixed 
number of turns of the bead ALk =15. This dynamics is 

performed with the global thermostat after thermal equilibrium 
has been reached. Note the diffusion of the plectoneme all along 
the DNA molecule. The DNA molecule is compos(;d of N = 300 
rigid cylinders of radius f = 2nm and / = 3.34nm. The other 
parameters of the simulation are given in Table 1 . This video is 
also available at http://vimeo.com/51918378. 
(AVI) 

Video S3 Simulated dynamics of a DNA molecule manipulated 

by magnetic tweezers with a stretching force / = 0.74 pN and fixed 
twisting torque r= 15 pN-nm. Global thermostat is used all over 
the simulation. At initial time the number of turns of the bead is 
zero and the torque starts to be applied to the bead (the rotation of 
the bead can be followed thanks to the small orange tag on the 
sphere). Because the applied torque (r= 15 pN-nm) is well above 
the buckling torque (r^ ~ 10 pN-nm), the DNA molecule is in the 
pure plectonemic state at the end of the simulation. The DNA 
molecule is composed of A = 300 rigid cylinders of radius r = 2 nm 
and / = 3.34 nm. The other parameters of the simulation are given 
in Table 1. This video is also available at http://vimeo.com/ 
51918151. 
(AVI) 
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